In this document we are going to calculate Evasion from the Vaisala Sensors and compare them to the eosFD flux measurements as well as see how they relate to other variables such as precipitation and temperature. First we import the data and for now we will filter so we are looking at non-injections
These new variables will help us plot each seperate day on a timescale of midnight -> midnight
df$day <- as.factor(substr(as.character(df$DateTime),1,10))
times <- substr(as.character(df$DateTime),12,16)
df$time <- strftime(df$DateTime, format="%H:%M")
df$time <- as.POSIXct(df$time, format="%H:%M")
Let’s visualize the changes in diel CO2 concentrations between a pair of sensors. Here we look at the entire reach, so stations 1 -> 4. We have a filter to remove some data that had issues.
filt <- df%>%
filter(!is.na(V1)&!is.na(V4))%>%
filter(DateTime > as.POSIXct("2019-07-13 00:00:00"))%>%
filter(DateTime < as.POSIXct("2019-08-13 19:00:00")|DateTime > as.POSIXct("2019-08-14 23:59"))%>%
filter(DateTime < as.POSIXct("2019-07-14 00:00")|DateTime > as.POSIXct("2019-07-14 23:59"))%>%
filter(!DateTime == as.POSIXct("2019-08-09 09:57:00"))
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-3200,-3200,-100,-100))
my_breaks <- c(0,15,50,166) # Create breaks
plot14Ppm <- ggplot(filt)+
geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
geom_point(aes(x = time, y = V4-V1,group = day, col=stn1_Q),lwd=2)+
scale_color_gradient(name = "Discharge (L/s)", trans = "log",
breaks = my_breaks, labels = my_breaks)+
labs(x = "Time",y="DCO2", title = "DCO2 Change station 4 minus station 1")+
theme_bw()
plot14Ppm
So now that we know what the change in CO2 between stations, it’s time to account for stream metabolism. Metabolism can be determined by taking the change in CO2 between stations before sunrise (photosynthesis is zero) and subtracting it from the change in CO2 at all other times. Let’s try to confirm this by looking the changes in DO over the same period.
Photosynthesis: Use CO2 and create O2 (1:1 ratio).
Respiration: Use O2, create CO2 (1:1 ratio).
Respiration can occur in the abscence of sunlight while photosynthesis cannot. At night, stream processes are respiration driven. During the day both respiration and photosynthesis occur. Is respiration constant throughout the day? if it is then pre-sunlight change in DO can be sen as respiration if no evasion of O2 is occuring. Therefore the area under the curve of the change in CO2 would be photosynthesis - respiration.
Metabolism
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-1,-1,1,1))
plot12DO <- ggplot(filt)+
geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
geom_point(aes(x = time, y = DO2_mg.L-DO1_mg.L,group = day, col=tempC_436),lwd=2)+
scale_color_gradient(name = "Discharge (L/s)", trans = "log",
breaks = my_breaks, labels = my_breaks)+
labs(x = "Time",y="DO (mg/L)", title = "Station 2 DO minus Station 1 DO")+
theme_bw()
plot12DO
## Warning: Removed 3136 rows containing missing values (geom_point).
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(0,0,4.5,4.5))
plot14DO <- ggplot(filt)+
geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
geom_point(aes(x = time, y = DO4_mg.L-DO1_mg.L,group = day, col=stn1_Q),lwd=2)+
scale_color_gradient(name = "Discharge (L/s)", trans = "log",
breaks = my_breaks, labels = my_breaks)+
labs(x = "Time",y="DO (mg/L)", title = "Station 4 DO minus Station 1 DO")+
theme_bw()
plot14DO
## Warning: Removed 3136 rows containing missing values (geom_point).
filt%>%
plot_ly(x = ~V1, y = ~DO1_mg.L,hoverinfo = "text",
text = ~paste("Time: ", DateTime,
"<br> Discharge: ", round(stn1_Q,2), "liters per second"))%>%
add_markers(color = ~tempC_436)%>%
layout(
title = "Station 1 Dissolved Oxygen vs CO2",
xaxis = list(title = "CO2 (ppm)", type = "log"),
yaxis = list(title = "Dissolved Oxygen (mg/L)")
)
## Warning: Ignoring 3136 observations
## Warning: Ignoring 3345 observations
## Warning: Removed 3136 rows containing missing values (geom_point).
## Warning: Removed 3136 rows containing missing values (geom_point).
There is a diel drop in DO in the late afternoon that seems to be related to temperature
## Warning: Removed 3136 rows containing missing values (geom_point).
(Station 4)
## Warning: Ignoring 3345 observations
Find the max change in pCO2 which will be when respiration is the least (before sunrise). This level serves as the baseline for that day (zero uptake). For every other time, we take the value of V1-V4 and subtract the zero uptake value for that day from it to determine uptake.
filtsub <- filt%>%
filter(time < as.POSIXct("06:30", format = "%H:%M"))%>%
group_by(day)%>%
mutate(baseDeltCo2 = max(V1-V4))%>%
ungroup()%>%
select(day,baseDeltCo2)%>%
distinct()
filt <- left_join(filt,filtsub)%>%
mutate(Uptake = baseDeltCo2-(V1-V4))
## Joining, by = "day"
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-500,-500,1800,1800))
plotUptake <- ggplot(filt)+
geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
geom_point(aes(x = time, y = Uptake,group = day, col=stn3_Q),lwd=2)+
scale_color_gradient(name = "Discharge (L/s)", trans = "log",
breaks = my_breaks, labels = my_breaks)+
labs(x = "Time",y="Uptake", title = "CO2 Respiration station 1 -> 4")+
theme_bw()
plotUptake
## Warning: Removed 78 rows containing missing values (geom_point).
Now that we know the change in CO2 and uptake, we can calculate evasion.
# Distances are 17,30,80
#filt$evasion_14 <- ((filt$V1-(filt$V4+filt$Uptake))/ 100) *.0018*(1000000/44.01)*filt$stn3_Q/1000 #with 'uptake'
filt$evasion_14 <- ((filt$V1-filt$V4)/ 100) *.0018*(1000000/44.01)*filt$stn3_Q/1000 #without
ggplot(filt, warning = FALSE, message=FALSE)+
geom_histogram(aes(x = evasion_14),fill = "#767676",color="#4B9CD3",binwidth = 1)+
geom_vline(xintercept = 0, linetype="dotted", color = "#4B9CD3", size=1.5)+
xlim(-10,40)+
ylim(0,2600)+
labs(x = "Evasion", y="Count of values",title = "Station 1 -> 4")+
theme_bw()
## Warning: Removed 2 rows containing missing values (geom_bar).
pptBar <- ggplot()+
geom_bar(data = pptDaily,stat = "identity", aes(x = DateTime, y = ppt_mm))+
geom_point(data = filt,aes(x = DateTime, y = evasion_14))
pptBar
##Evasion
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-35,-35,40,40))
plot14 <- ggplot(filt)+
geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
geom_line(aes(x = time, y = evasion_14,group = day, col=stn1_Q),lwd=2)+
scale_color_gradient(name = "Discharge (L/s)", trans = "log",
breaks = my_breaks, labels = my_breaks)+
labs(x = "Time",y="Evasion", title = "Evasion station 1 -> 4")+
theme_bw()
plot14
ggplot(filt)+
geom_point(shape=21,aes(x = evasion_14, y=Flux_1, fill = "eosFD 1"), fill="#4B9CD3", alpha =.6,size=2)+
geom_point(shape=21,aes(x = evasion_14, y=Flux_2, fill = "eosFD 2"), fill="#C7DA7D", alpha = .6,size=2)+
geom_smooth(aes(x=evasion_14,y=Flux_1),method='lm',formula=y~x, color = "#4B9CD3")+
geom_smooth(aes(x=evasion_14,y=Flux_2),method='lm',formula=y~x,color ="#C7DA7D")+
xlim(0,36)+
ylim(0,6)
## Warning: Removed 3394 rows containing non-finite values (stat_smooth).
## Warning: Removed 3379 rows containing non-finite values (stat_smooth).
## Warning: Removed 3394 rows containing missing values (geom_point).
## Warning: Removed 3379 rows containing missing values (geom_point).
Using the USGS streammetabolizer package, daily estimates of GPP and ER can be calculated at each of the DO sensor locations. This is done in a seperate script located in: ‘Ecuador/Analysis/Stream_Metabolism/metabolism.R’. The results are used to estimate daily evasion for the reach. The values expressed in the results of the streammetabolizer tool are grams of O2 / m^2 /d^-1. We need to convert this into uMols and then multiply by -1 to determine uMols CO2.
pred1 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_1_Predictions.csv"))%>%
select(date,GPP,ER)%>%
mutate(Station = "1")
pred2 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_2_Predictions.csv"))%>%
select(date,GPP,ER)%>%
mutate(Station = "2")
pred4 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_4_Predictions.csv"))%>%
select(date,GPP,ER)%>%
mutate(Station = "4")
metab <- rbind(pred1,pred2,pred4)%>%
mutate(DateTime = as.POSIXct(date))
metab%>%
plot_ly(x = ~DateTime, y = ~GPP,hoverinfo = "text",
text = ~paste("Station: ", Station,
"<br> Date: ", DateTime,
"<br> ER: ", round(ER,2)))%>%
add_markers(color = ~Station)%>%
layout(
title = "Gross Primary Productivity (DO)",
xaxis = list(title = "Date"),
yaxis = list(title = "GPP (g/m^2/d^-1)")
)
## Warning: Ignoring 42 observations
1 mole of CO2 = 31.9988 grams. We average the GPP and ER values from the three stations to estimate GPP and ER for the entire reach.
metabCO2 <- metab%>%
mutate(GPP_CO2uMol = GPP/31.9988*-1000000,
ER_CO2uMol = ER/31.9988*-1000000)%>%
group_by(DateTime)%>%
mutate(AvgGPP_CO2 = mean(GPP_CO2uMol),
AvgER_CO2 = mean(ER_CO2uMol))
metabCO2%>%
plot_ly(x = ~DateTime, y = ~GPP_CO2uMol,hoverinfo = "text",
text = ~paste("Station: ", Station,
"<br> Date: ", DateTime,
"<br> ER: ", round(ER_CO2uMol,2)))%>%
add_markers(color = ~Station)%>%
layout(
title = "Gross Primary Productivity (uMols CO2)",
xaxis = list(title = "Date"),
yaxis = list(title = "GPP (uMols/m^2/d^-1)")
)
## Warning: Ignoring 42 observations
metabCO2%>%
plot_ly(x = ~DateTime, y = ~ER_CO2uMol,hoverinfo = "text",
text = ~paste("Station: ", Station,
"<br> Date: ", DateTime,
"<br> ER: ", round(GPP_CO2uMol,2)))%>%
add_markers(color = ~Station)%>%
layout(
title = "Ecosystem Respiration (uMols CO2)",
xaxis = list(title = "Date"),
yaxis = list(title = "ER (uMols/m^2/d^-1)")
)
## Warning: Ignoring 42 observations
evasionDaily <- filt%>%
group_by(day)%>%
mutate(DailyFlux = mean(evasion_14*86400))%>%
select(day,DailyFlux)%>%
distinct()%>%
rename("date"="day")
metabAvg <- metabCO2%>%
select(date,AvgGPP_CO2,AvgER_CO2)%>%
distinct()
## Adding missing grouping variables: `DateTime`
adjustEvasion <- metabAvg%>%
left_join(evasionDaily,by="date")%>%
mutate(AdjEvasion = DailyFlux+AvgGPP_CO2+AvgER_CO2)%>%
select(-date)
## Warning: Column `date` joining factors with different levels, coercing to
## character vector
adjustEvasion%>%
plot_ly(x = ~DateTime, y = ~AdjEvasion,hoverinfo = "text",
text = ~paste("Evasion: ",round(AdjEvasion,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>",
"<br> GPP: ", round(AvgGPP_CO2,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>",
"<br> ER: ", round(AvgER_CO2,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>"))%>%
add_markers(color = "Evasion", size=50)%>%
add_markers(y = ~AvgER_CO2, color = "Respiration", size = 50)%>%
add_markers(y = ~abs(AvgGPP_CO2), color = "Photosynthesis", size=50)%>%
layout(
title = "Daily Evasion (Entire Reach)",
xaxis = list(title = "Date"),
yaxis = list(title = "uMols m<sup>2</sup>d<sup>-1</sup>", type = "log")
)
## Warning: Ignoring 15 observations
## Warning: Ignoring 14 observations
## Warning: Ignoring 14 observations